Time reversal symmetry breaking superconductivity in the honeycomb t-J model 
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We report the theoretical discovery of a novel time reversal symmetry breaking superconducting 
state in the t-J model on the honeycomb lattice, based on a recently developed variational method 
- the Grassmann tensor product state approach. As a benchmark, we use exact diagonalization 
(ED) and density matrix renormalization (DMRG) methods to check our results on small clusters. 
Remarkably, we find systematical consistency for the ground state energy as well as other physical 
quantities, such as the staggered magnetization. At low doping, the superconductivity coexists with 
anti-ferromagnetic ordering. 
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Introduction 

Since the discovery of high-temperature superconduc- 
tivity in the cuprates[l], many strongly correlated mod- 
els have been intensively studied. One of the simplest of 
these is the t-J model[2]: 

H t -j = t A/ c 3,° + h - c - + J Y1 ( S i ■ S J - \^ n j) > 

(ij) V + J 

(1) 

where c\ a is the electron operator defined in the no- 
double-occupancy subspace. This model can be derived 
from the strong-coupling limit of the Hubbard model. 
Despite its simplicity and extensive study, the nature of 
the ground states of Eq. (1) is still controversial. 

A strong correlation view of the t-J model was ad- 
vanced by Anderson, who conjectured the relevance of a 
Resonating Valence Bond (RVB) state[3] as a low energy 
state for Eq. (1) when doped. When undoped, the RVB 
state is a spin singlet, with no symmetry breaking, and 
describes a "quantum spin liquid" . At low temperature 
the mobile carriers in the doped RVB state behave as 
bosons and condense, forming a state indistinguishable 
in terms of symmetry from a singlet BCS superconduc- 
tor. A further development was the introduction of a 
projected mean-field wavefunction - the projection re- 
moving all components of the wavefunction with doubly 
occupied sites which could be used variationally[4, 5] 

Presently, this variational method remains one of the 
few numerical tools to study t-J like models which work 
directly at T — and can treat significant system sizes. 
However, due to the special form of the variational wave- 
function, one may be concerned about bias: very general 
states in the low energy subspace cannot be investigated. 

Recently new numerical tools have been developed to 
investigate much more general low energy states in t- 
J like models beyond the projective method. One novel 



construction builds Tensor Product States (TPS's)[7-12], 
which can be conveniently studied and have been applied 
to many spin systems [8-14]. This new class of variational 
states do not assume any specifical ordering pattern a 
priori and can describe very general states as long as 
their entanglement entropy satisfy perimeter law. Re- 
cently, this method has been generalized to fermionic 
systems [15-20]. Among many different generalizations, 
the Grassmann Tensor Product States (GTPS's)[20] were 
shown to be closely related to projective states. They are 
able to describe a class of projective wavefunctions faith- 
fully, including in particular the short-range RVB states. 

The TPS and GTPS construction is particularly well 
suited to trivalent lattices, the simplest of which is the 
honeycomb lattice in two dimensions. Like the square lat- 
tice appropriate for the cuprates, the honeycomb lattice 
is bipartite and naturally supports an antiferromagnctic 
(AF) state at half-filling in the strong coupling (Heisen- 
berg) limit. Similarities to cuprate physics may be ex- 
pected. Moreover, several numerical studies have identi- 
fied a possible quantum spin liquid state on this lattice 
at half-filling when additional quantum fluctuations are 
included in the Hubbard[6] and Heisenberg[21, 22] mod- 
els. Thus the doped t-J model on the honeycomb lattice 
seems a promising venue to explore RVB ideas. 

Here, we investigate the ground state of this system us- 
ing the recently developed GTPS approach. The GTPS 
results are benchmarked by comparison with exact di- 
agonalization (ED) and density matrix renormalization 
group (DMRG) calculations. Our results are systemati- 
cally consistent with these non-variational, exact meth- 
ods on small clusters. The principle result of the GTPS 
calculations is that the ground state at non-zero doping is 
a time reversal symmetry breaking d + id wave supercon- 
ductor. Some physical rationale for this result are given 
at the end of this Letter. 



2 



9 




e„ de„ de„. e„, 



FIG. 1: Graphic representation of the GTPS on a honey- 
comb lattice. Ta and Tb, which contain 9 are defined on 
the sublattices A and B for each unit cell. The Grassmann 
metric g containing dO is defined on the links that connect the 
Grassmann tensors and Tg . The blue lines represent the 
fermion parity even indices while the red lines represent the 
fermion parity odd indices of the virtual states. Notice an ar- 
row from A to B represents the ordering convention d6 a d9 a i 
that we use for the Grassmann metric. 
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Results 

The variational ansatz: We use the standard form 
of GTPS as our variational wavefunction. We further as- 
sume a translationally invariant ansatz, and thus is spec- 
ified by just two different Grassmann tensors Ta,Tb on 
sublattice A, B of each unit cell: 
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We notice the symbol tTr means tensor contraction of the 
inner indices {a}. Here 9 a {fl,i)i^a{fl,i) are the Grass- 
mann numbers and dual Grassmann numbers, respec- 
tively, defined on the link a(b, c), which satisfy the Grass- 
mann algebra: 



V a V0 — —UpUa, 

d9 a 6p = 5 a p 



dO a d9 p = -d6pd6 a , 
d9 a l = 0. (4) 



As shown in Fig.l, a, b, c — 1,2,...,!? are the virtual 
indices carrying a fermion parity P*{a) =0,1. In this 
paper, we choose D to be even and assume there are equal 
numbers of fermion parity even/odd indices, which might 
be not necessary in general. Those indices with odd par- 
ity are always associated with a Grassmann number on 
the corresponding link and the metric g aa / is the Grass- 
mann generalization of the canonical delta function. The 
complex coefficients T n A 
tional parameters. 
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FIG. 2: Ground state energy as a function of doping. As a 
benchmark, we performed the ED calculation and DMRG cal- 
culations for small system size. Insert: Stagger magnetization 
as a function of doping. 



Notice that is the physical index of the t-J model 
on site i, which can take three different values, o, t and 
I, representing the hole, spin-up electron and spin-down 
electron states. We choose the hole representation in our 
calculations, and thus the hole state has an odd parity 
Pf (o) = 1 while the electron states have even parity pf (f 
,4-) = 0. On each site, the non-zero components of the 
Grassmann tensors should satisfy the parity conservation 
constraint: 

P f (m,) + P f (a) + P f (6) + P f (c) = 0(mod 2). (5) 

Since the wavefunction Eq.(3) does not have a definite 
fermion number, we use the grand canonical ensemble, 
adding a chemical potential term to Eq.(l) to control 
the average hole concentration. 

We use the imaginary time evolution method[23] to 
update the GTPS from a random state. Then we use 
the weighted Grassmann-tensor-entanglcment renormal- 
ization group (wGTERG) mcthod[20, 23] to calculate 
physical quantities. (See methods section for details.) The 
total system size ranges up to 2 x 27 2 sites and all calcula- 
tions are performed with periodic boundary conditions. 
The largest virtual dimension of the GTPS considered 
was 12. To ensure convergence of the wTERG method, 
we kept D cut (defined in Refs. [20, 23]) up to 130 for 
D = 4, 6, 8, 10, which gave relative errors for physical 
quantities of order 10 -3 . For D = 12, we kept D cut up 
to 152. 

Ground state energy and staggered magneti- 
zation: At half-filling, the t-J model reduces to the 
Heisenberg model. In this case, we find the converged 
ground state energies per site are —0.5439 for D = 10 
and -0.5441 for D = 12(the term 
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here), which is consistent with previous TPS study[12] 
(with virtual dimension D — 5 and 6, since all the com- 
ponents of GTPS with odd fermion parity virtual indices 
vanish in this case) and a recent Quantum Monte Carlo 
(QMC) result E = -0.54455(20). Despite good agree- 
ment with the ground state energy, the staggered mag- 
netization m obtained from our calculations is larger than 
the QMC result m = 0.2681(8). We hnd m = 0.3257 for 
D = 10 and m = 0.3239 for D = 12, also consistent with 
the previous TPS study[10-12]. 

Nevertheless, we emphasize that the variational ap- 
proach indeed obtains the correct phase. Actually, a re- 
cent study for square lattice Heisenberg model shows m 
can be consistent with the Quantum Monte Carlo (QMC) 
result if D is sufficient large [24] . 

Much more interesting physics arises after we dope the 
system. (We consider t/J = 3). As seen in Fig. 2, the 
ground state energy shows a marked increase in D de- 
pendence with increasing hole doping S. As a bench- 
mark, we perform the ED and DMRG calculations for 
small periodic clusters with N sites (N = 18 for ED and 
N = 36, 54 for DMRG). These two methods are the only 
unbiased methods for frustrated systems that avoid the 
sign problem, but are restricted to relatively small sys- 
tems. To ensure the convergence of the DMRG, we keep 
up to 8000 states and make the truncation errors less 
than 10~ 9 in our N = 54 calculations. Up to D = 10, 
we find a systematic convergence of the ground energy. 
Some data points around S = 0.1 for D = 12 have slightly 
higher energy than D — 10, because D cut = 152 is still 
not large enough for convergence at D = 12. 

As shown in the insert of Fig. 2, the staggered magneti- 
zation to has an even larger D dependence than energy at 
finite doping. However, up to D = 12, the data appears 
to converge to a relatively well-defined curve indicating 
vanishing AF order for 6 > 0.1. 

Superconductivity: Next we turn to the interesting 
question of whether the doped antiferromagnetic Mott 
insulator on the honeycomb lattice supports supercon- 
ductivity, and if so, what its pairing symmetry is. To 
answer this, we calculate the real space superconduct- 
ing (SC) order parameters in the spin singlet channel 
A s = (ci^Cj^ — Ci.^Cj^), where i and j are nearest 
neighboring sites. Because we use a chemical potential to 
control the hole concentration, the charge U(l) symmetry 
can be spontaneous broken in the variational approach, 
which allows A s to be directly measured rather than 
through its two-point correlation function. As shown in 
the main panel of Fig. 3, up to S = 0.15 we find a non- 
zero singlet SC order parameter for the whole region. 
Strikingly, we find that the SC state breaks time reversal 
symmetry. By measuring the SC order parameters for 
the three inequivalent nearest-neighbor bonds, we found 
A s a /A s b = A S JA S C = A s JA s a = e 10 with 9 = ±^f. This 
pairing symmetry is usually called d + id wave. To ex- 
clude the possibility that this results from trapping in 
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FIG. 3: SC order parameters as a function of doping. 



a unstable local minimum, we repeated the calculations 
with many different random tensors, and all cases con- 
verged to the same results. Moreover, we also checked 
that the SC order parameter vanishes at large t/J to 
make sure that the existing of SC order is consequence of 
spontaneous symmetry breaking. (For D — 6, the critical 
value is around 15 at 5 ~ 0.3 ) 

The existence of SC order is observed in our numerical 
study up to 5 — 0.4. However, a much larger inner di- 
mension D is required for the convergence of ground state 
energy at larger hole concentration, which is beyond the 
scope of this paper. (The GTPS variational ansatz we use 
in this paper is designed to help us understand the na- 
ture of Mott physics; at large doping, the Mott physics 
becomes less important and can be studied much better 
by other methods.) 

Coexisting phases at low doping: Interestingly, 
we find the SC and AF order coexisting in the regime 
< S < 0.1. A physical consequence of the microscopic 
coexistance is that triplet pairing is induced. The in- 
sert of Fig. 3 shows the amplitude of the triplet order 
parameter as a function of doping. Since the triplet pair- 
ing order parameter has three independent components 
A t — -^Ci^ a (ia v <j) a f3Cj y p = de I(?i (Here i G A, j e B and 
<p is the phase of SC order parameter.), we can define the 



amplitude of triplet order parameter as A t = y A£ • A t . 

The phase shift of 27r/3 on the three inequivalent bonds 
is also observed for all the triplet components. We fur- 
ther check the internal spin direction of the triplet d 
vector and find it is always anti-parallel to the Neel 
vector (SnccI = (Si) — (Sj)). At larger doping 5 > 0.1, 
the triplet order parameter has a very strong D depen- 
dence, so at present we are unable to determine whether 
it ultimately vanishes or remains non-zero in the D — > oo 
limit. We leave this issue for future work. By fully using 
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FIG. 4: The comparison of DMRG/ED and GTPS calcula- 
tions for stagger magnetization m. 

all symmetry quantum numbers and other techniques like 
high performance simulation on GPUs, we can in princi- 
ple handle D up to 20 — 30. 

Discussion 

The calculations in this paper by construction ig- 
nore any possible phases breaking translational sym- 
metry. Natural candidates are spiral or striped anti- 
ferromagnetic phases, as have been heavily discussed for 
the square lattice. A weak coupling perspective suggests 
this may be unlikely. 

In the weak-coupling limit of Hubbard model on hon- 
eycomb lattice, the system is a semimetal with two Dirac 
cones. With a non-zero Hubbard U, commensurate AF 
fluctuations at zero momentum manifest in the inter- 
band susceptibility, becoming stronger and stronger with 
increasing U. At sufficiently large U commensurate AF 
order develops (recent numerics find a narrow region of 
intermediate spin liquid phase [6]). At small but finite 
doping, the Dirac cones become pockets. In this case, 
the total intra-band spin susceptibility shows a constant 
behavior for small q(q < 2fc/)[25] while the total inter- 
band spin susceptibility still peaks at zero momentum. 
Thus, we argue the commensurate AF fluctuations at 
zero momentum still dominate for sufficient small hole 
concentration. 

Our DMRG and ED calculations also support this ar- 
gument. Up to 54 sites, we do not find any evidence of 
incommensurate spin-spin correlation. Fig. 4 shows the 
comparison of DMRG calculations and GTPS calcula- 
tions for staggered magnetization m on small systems. 
We find the GTPS results are comparable with extrapo- 
lations of m for infinite size systems (the GTPS results 
for m are somewhat larger close to half-filling due to in- 
sufficient tensor dimension D). This agreement also im- 



plies our variational wavefunction is very close to the true 
ground state. 

A more general concern is whether the GTPS tends 
to overestimate SC order at large doping, due to its 
non-conservation of charge(note that the projected wave- 
function approach has a rather strong tendency to pro- 
duce superconducting states). The observed SC order 
parameter is "small" in terms of the natural upper limit 
(cc) < 0.1(5 6. Nevertheless, our results are consistent 
with other approaches: (a) in the mean field theory for 
the honeycomb t — J model, only the d + id pairing chan- 
nel gains energy [26], and (b) in the weak coupling limit 
of Hubbard model, very recent renormalization group 
studies also find d + id superconductivity around quarter 
filling[27-29]. 

In conclusion, we report the theoretical discovery of 
a d + id wave superconducting state in t-J model on a 
honeycomb lattice, based on a recently developed vari- 
ational method - the GTPS approach. At low doping, 
AF ordering coexists with the SC ordering. In the co- 
existance regime, a spin triplet pairing with the same 
phase shift is induced and its triplet d vector is anti- 
parallel with the Neel vector. It would be interesting 
to search for this physics in experiment. The recently 
discovered spin 1/2 honeycomb lattice antiferromagnet 
InV!/3Cu2/3O3[30] would be an appealing candidate if it 
could be doped. 

Methods 

The imaginary time evolution algorithm of 
GTPS: In this paper, we use the (simplified) imaginary 
time evolution method[23] to update the GTPS varia- 
tional wave function. In the hole representation, we can 
decompose the c\ as Cj !<T = ft|&i ;o -. Here the holon hi 
is a fermion while the spinon b^ a is a boson. The no- 
double-occupancy constraint reads: 

E 6 tA- + ^ = 1 ( 6 ) 



Under this representation, we can rewrite the t — J model 
as: 




(7) 

where 

Due to the no-double-occupancy constraint Eq. (6), the 
spin up/down states | t = &l-|~n\|0) and the hole 

state |o) = h\\0) form a complete basis for each site. The 
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FIG. 5: A schematic plot for the imaginary time evolution 
algorithm. 



closure condition reads: 



ItiXti \ + \U)(U \ + \oi)(oi\ = 1 



(9) 



As already has been discussed in Ref. [23], we need to 
use the fermion coherent state representation to perform 
the imaginary time evolution algorithm for GTPS. Let 
us introduce the fermion coherent state of holon \rji) = 
|0) — r7j/i||0).(77j is a Grassmann variable here.) In this 
new basis, the closure relation Eq. (10) becomes: 

I tiXti \ + \U) (U I + / drhdrhh) (m 1 = 1 (10) 

The variational ground state can be determined through 
imaginary time evolution: 



tf G } = e- rH *-'|*o), t^oo 



(11) 



For sufficient thin time slice St, we can decompose 
e -6rH t -j ag . 

g-Jrflw „ e -«Tfl?_ Je -«rfl*_ Je -«rff/_ J (12) 

Here x, y, z represent three different directions of the hon- 
eycomb lattice and each H^j^ term contains a summa- 



t-J 



tion of non-overlapped two body Hamiltonian H\ 

hij • Thus, for sufficient thin time slice, we can 
apply the evolution operator along x(y, z) direction sep- 
arately. 

In the fermion coherent state representation, the two 
body evolution operator along x direction e~ Shi i can be 
expressed as a Grassmann valued matrix acting on the 
two sites Hilbert space: 



--E' 



^{n j ) pS{mi \r li ) pl ^\n[f { <\i j f t ^ (13) 



Here rrtj is the fermion coherent states indices valued as 
t, 4- or V while vrii is the usual Fock state indices valued 
as f, I or o. We notice the spin states remain the same 
in both representations. 

By applying the method developed in Ref. [23] , we can 
successfully update the (complex) variational parameters 



TTL and T" 3 ,,, , as: 

A:abc B:a b c 



rplrrii _ ( \P f (mi)P f (a) 

1 A.„U„ ~ I, ) 



A:abc 




bcrrii 



T'" lj - ( \P , (m i )PUa') V J V TZ. , , 

V A b>V A c> 



where U and V are determined by the singular value 
dccomposition(SVD) of the following matrix M: 



M 



bcjni :b' c'jn 



x (_)[ p/ ( m 'i) + pf i m 'j)]P S i a )^_^ pf ( m 'j)[ pS ( m 'i) + pf {b)+P f (c)] 
X \ ) >£j m ' i m' :j 1 A;abc 1 B;ab'c' \ L0 ) 

We keep the largest Dth singular values: 

D 

^^-bcrrii\b ! c ! nij — ^ ^ Ubcrrij \a, A aVb' c'vrij -a (1^*) 



Similar as in the usual TPS case[10], the environment 
weight vectors K x ^ y - Z ^ can be initialized as 1 and then 
updated during the time evolution. For example, A x is 
updated as A' in the above evolution scheme. 

The GTERG / wGTERG algorithm: After deter- 
mining the variational ground state of GTPS by per- 
forming the imaginary time evolution, we can compute 
the physical quantities by using the GTERG / wGTERG 
method developed in Ref [20, 23]. This method can be 
regarded as the Grassmann variable generalization of 
the usual TERG/wTERG method [9, 23], where a coarse 
graining procedure is designed to calculate the physical 
quantities of TPS efficiently. In Fig. 7, we plot the vari- 
ational ground state energy(D — 10) as a function of 
doping with different D cut (Number of eigenvalues kept 
in the wGTERG algorithm). We find the relative error 
is of order 10 -3 . 
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FIG. 6: A schematic plot for the renormalization algorithm on 
honeycomb lattice. Similar as the simplified imaginary time 
evolution algorithm, we use a weighting vector V" to mimic 
the environment effect. As has been discussed in Ref.[23], the 
initial value of Q can be determined by A on the corresponding 
link and will be updated during the RG scheme. 
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FIG. 7: Variational ground state energy(D = 10) as a func- 
tion of doping with different D cu t (Number of eigenvalues kept 
in the wGTERG algorithm). It is shown all the data points 
almost collapse on the same curve. We find the relative error 
is of order 10 -3 . 
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